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ABSTRACT 


The  optimal  minimum  time  control  (i.e.  bang-bang  controller)  is  applied  to 
the  fast  reaction  missile  defense  problem.  From  Pontryagin,  the  optimal  control 
was  determined  to  be  a  function  of  the  adjoint  in  the  minimization  of  the 
Hamiltonian.  The  control  may  also  be  posed  either  as  a  function  of  time  or  as  a 
function  of  the  states.  The  state  space  can  be  partitioned  into  regions,  surfaces 
and  curves  where  the  optimal  control  action  is  either  its  maximum  plus  or  minus 
N. 

In  missile  simulation  problems,  the  method  of  adjoints  is  often  used  in 
parametric  studies  of  errors  and  miss  distance.  This  technique  is  developed 
both  graphically  and  mathematically,  and  is  used  here  to  help  one  visualize  the 
solution  trajectory  and  families  of  optimal  trajectories  for  all  possible  initial 
conditions. 
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I.  INTRODUCTION 


In  this  era  of  technological  advancement,  more  and  more  effort  is  being 
devoted  to  automation  to  increase  speed  and  efficiency.  Everything  is 
becoming  faster,  smaller,  more  efficient  and  more  effective,  from  fuzzy  logic 
controlled  appliances  to  smart  weapon  systems.  As  our  engineering  and  design 
of  physical  systems  improves  we  are  able  to  generate  faster  and  more  accurate 
mechanical  devices,  and  the  control  systems  for  such  devices  must  improve  as 
well.  In  pushing  »o  develop  the  fastest,  most  efficient  controller  foi  whatever 
our  application,  we  introduce  the  bang-bang  controller. 

Many  of  our  current  defensive  missile  systems  were  designed  to  be  used 
against  threats  that  are  no  longer  of  greatest  importance.  A  point  defense 
missile  may  now  be  required  to  bring  down  a  target  that  has  a  significant  speed 
advantage,  and  be  able  to  do  it  in  such  a  way  as  to  control  the  fallout  after  the 
impact. 

Even  if  we  redesign  our  systems,  current  economic  conditions  make  it 
unlikely  that  we  could  build  up  an  arsenal  of  new  weapon  systems.  However, 
we  could  upgrade  our  currently  existing  arsenal  by  replacing  the  original  control 
logic  with  a  newer,  more  capable  controller. 

In  this  report  we  will  develop  second  and  third  order  minimum  time  optimal 
controllers  and  apply  them  to  a  fast  reaction  missile  defense  problem.  While  we 
are  using  the  missile  defense  problem  as  our  example,  it  shoul*^.  be  noted  that 
the  controllers  are  not  limited  to  this  example  and  may  be  applied  as  generic 
controllers  for  a  wide  variety  of  systems. 
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II.  SECOND  ORDER  CONTROLLER 


A.  MINIMUM  TIME  OPTIMAL  CONTROL 

Designing  a  control  system  is  frequently  a  hit-and-miss  process  using  a 
variety  of  design  techniques  to  iteratively  create  a  system  that  meets  specific 
criteria.  The  performance  of  the  system  may  be  defined  in  the  time  and 
frequency  domain  in  terms  of  overshoot,  settling  time,  etc.,  or  it  may  be  defined 
by  some  externally  measured  criteria.  For  example,  a  satellite  control  system 
may  be  designed  to  minimize  fuel  consumption. 

"The  objective  of  optimal  control  theory  is  to  determine  the  control  signals 
that  will  cause  a  process  to  satisfy  the  physical  constraints  and  at  the  same 
time  minimize  (or  maximize)  some  performance  criterion."  jl] 

1.  Problem  Definition 

The  system  and  optimization  problem  is  defined  as 

x  =  Ax  +  Bu  (2.1) 

with  x(0)  =  c,  and  x(tf)  =  0,  while  minimizing 


‘f 

J  =  Jdt.  (2.2) 

0 

Let  us  consider  the  simple  example  with 


'0 

r 

"O' 

X  = 

0 

0 

x  + 

1 

(2.3) 


From  Pontryagin  [1 1,  we  find  we  can  minimize  J  by  minimizing  the  Hamiltonian 


H  =  1  +  p,X2  +  P2U. 


(2.4) 
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This  is  minimum  when  u  is  operating  at  its  maximum  possible  value  and  with 
opposite  the  sign  of  the  adjoint  p2.  Thus  we  have 

u  =  -Nsign(p2)  (2.5) 


where 


This  has  a  solution 


P 


^=r  ^ 

<3x  -1 


Pl=Ci 

P2  =-Cit  +  C2. 

A  typical  solution  would  appear  as  seen  in  Fig.  2. 1 . 


(2.6) 


(2.7) 

(2.8) 


Figure  2.1  Solution  to  Minimized  Hamiltonian 

Note  that  from  the  adjoint  solution,  the  control  may  change  sign  only 
once.  We  would  like  to  solve  this  problem  for  all  possible  initial  conditions  and 
only  one  terminal  condition  (x(tf)  =  0).  Hence  it  makes  sense  to  look  at  this 
problem  in  negative  time  (adjoint),  starting  from  the  end  point  of  motion.  From 
the  uniqueness  theorem  in  ordinary  differential  equations,  only  two  trajectories 
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can  emanate  from  the  origin;  one  for  u  =  -N  and  one  for  u  =  +N,  as  shown  in 
Figure  2.2.  In  this  second  order  example,  our  solution  is  constrained  to  the  xi, 
X2  plane.  These  negative  time  trajectories  divide  the  state  space  (here  a  plane) 
into  two  parts.  Solution  trajectories  emanating  from  these  curves  constitute  all 
possible  solutions  for  all  possible  initial  conditions. 


Figure  2.2  Zero  Trajectory  Curves 
2.  Minimum  Time  Trajectories,  Second  Order 
Consider  the  system 


X  = 


0  1 
0  0 


x  + 


u. 


This  can  be  represented  in  flow  diagram  form  as 


(2.9) 
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The  control,  u,  may  be  defined  for  a  bang-bang  control  system  as  ±N,  where  N 
is  some  constant.  Given  any  starting  values  for  the  states,  there  are  only  two 
possible  paths  for  the  states  to  take;  one  corresponding  to  u  =  +N,  and  one 
corresponding  to  u  =  -N.  If  we  desire  to  drive  all  states  to  zero,  and  we  know 
that  we  can  only  apply  u  =  ±N,  Figure  2.3  shows  the  paths  followed  by  the 
states  given  various  starting  values. 


Figure  2.3  Minimum  Time  Trajectory  Solution  Curves 
3.  Second  Order  Solution  Trajectories 
This  system  has  a  solution 

x(t)  =  0(t,O)x(O)-i-  A(t,0)u(0)  (2.10) 

where 

(2.11) 

and 
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(2.12) 


A= r 

Jo 


Expanding  ,  we  get 
e^‘  = 


1  + At  +  - 

1  A  2  2 

—  At 

L 

2! 

'1  0‘ 

0  r 

0  1 

+ 

0  0 

1  O' 

ro  ti 

““ 

0  1 

+ 

1 

o 

o 

1 

t  +  — 
2! 

1 

+ - 

2! 


'0 

1]' 

_0 

oJ_ 

0  1 
0  0 


— 1 
o 

o 

'1  f 

t  +...= 

o 

o 
_ 1 

0  1 

(2.13) 


Since  is  the  zero  matrix,  all  higher  orders  of  A  will  also  be  the  zero  matrix 
and  we  may  drop  them.  Evaluating  A  (for  u  fixed)  we  find 


A  =re^(‘-^)Bdx=r 
Jo  Jo 


1 — 

1 

1 _ 

'O' 

Lo  1  J 

1 

dx 


=i: 

The  solution  for  this  system  is 


't-x' 

dx  = 

J 

1 

1-  -I 

1 

= 

2 

X 

0 

t 

(2.14) 


x(t)  = 


'1  t 
0  1 


x(0)+ 


■j-t^ 
2  '■ 


U(0), 


(2.15) 


or  in  scalar  equations 

Xi(t)  =  Xi(0)+X2(0)t  +  iu(0)t^  (2.16) 

X2(t)  =  X2(0)  +  u(0)t.  (2.17) 

These  equations  describe  the  states  as  a  function  of  time  given  any  initial 
conditions  and  the  fixed  control  effort,  u. 


B.  ANALYTIC  SOLUTION  FOR  SECOND  ORDER  SWITCHING 
TIME 

We  may  now  treat  this  system  as  a  boundary  value  problem  and 
analytically  solve  for  switching  times  of  the  control  effort. 
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1.  Solving  Boundary  Value  Problems 

Since  the  control  is  piecewise  constant,  (±N),  we  can  separate  the 
problem  into  two  pieces  and  match  boundary  values  at  the  point  where  u 
changes  sign.  In  other  words,  if  the  system  moves  from  point  x(to)  through 
point  x(ts)  to  point  x(tf),  we  can  solve  the  problem  in  two  parts.  Each  of  our 
boundary  value  problems  can  be  stated  in  such  a  way  as  to  supply  simplifying 
boundary  conditions,  i.e.,  setting  initial  or  final  values  to  zero.  Optimal  control 
theory  tells  us  that  for  a  second  order  system  there  will  be  at  most  one 
switching  time,  the  change  from  u  =  +N  to  u  =  -N,  or  visa  versa.  Let  us 
consider  first  the  solution  for  our  system  with  some  initial  condition  [xi(0), 
X2(0)],  as  shown  in  Figure  2.4. 


Figure  2.4  Minimum  Time  Trajectory  From  a  Fixed  Initial 
Condition 

The  time  from  t=0  to  t*  is  the  length  of  time  before  the  control  effort 
changes  sign.  We  cannot  solve  directly  for  t*  because  we  do  not  have  enough 
information  on  the  boundary  values.  Since  we  are  solving  this  problem  for 
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arbitrary  initial  conditions,  Xi(0)  and  X2(0),  xi(ts)  and  X2(ts)  are  unknown.  We 
do  know,  however,  that  the  final  states,  Xi(tf)  and  X2(tf),  will  be  zero,  and  from 
this  we  can  determine  the  zero-trajectory  curve  in  negative  time  from  the 
origin.  This  curve,  shown  in  Figure  2.5.,  is  the  only  curve  that  will  go  through 
the  origin  with  u  =  +N.  Therefore  the  switching  time  will  occur  when  the  path 
of  a  u  =  -N  curve  intersects  this  curve.  There  are  an  infinite  number  of  u  =  -N 
curves  intersecting  this  curve,  however,  only  one  curve  will  go  through  any 
particular  set  of  initial  conditions. 

To  simplify  the  problem,  we  will  start  with  initial  conditions  on  the  xi 
axis.  Define  Xi(0)  =  Ci,  X2(0)  =  0,  xi(tf)  =  X2(tf)  =  0,  and  u  =  -N.  This  situation 
is  described  in  Figure  2.6  where  ti  may  be  defined  as  the  time  from  to  to  t*. 
Equation  (2.16)  may  now  be  written 

Xj(t)  =  Xj(0)  +  X2(0)t  +  ^ut^  =  Xi(0)-^Nt^  (2.18) 


Since  the  curves  for  u  =  +N  and  u  =  -N  are  symmetric,  it  can  be  noted  that 


Xi(ti)  =  '5X,(0) 

(2.19) 

and  therefore 

ix,(0)  =  x,(0)-iNt,* 

(2.20) 

and  finally 

,  _  ^i(O) 

(2.21) 

This  is  valid  for  any  initial  conditions  such  that  xi(0)  >  0,  X2(0)  =  0  and  u 

=  -N. 
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Figure  2.5  Simplified  Boundary  Value  Problem  #1 
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The  next  simplification  involves  getting  from  the  positive  X2  axis  down 
to  the  positive  xi  axis  as  shown  in  Figure  2.7.  Given  xi(0)  =  0,  X2(0)  =  positive 
real,  Xi(tf)  =  positive  real,  X2(tf)  =  0,  and  u  =  -N. 

Equation  (2.17)  may  be  written 

X2(t)  =  X2(0)-Nt.  (2.22) 


At  t  =  t2  ,  X2(t2)  =  0  giving 


0  =  x,(0)-Nt2 


(2.23) 


,  _X2(0) 


(2.24) 


The  final  stage  is  to  translate  the  initial  condition  off  of  the  X2  axis  to 
some  unknown  initial  condition.  The  time  required  for  the  states  to  get  from 
some  initial  condition  X2(0)  to  X2(t)  is  not  based  on  the  xi  state,  and  is  therefore 
always  equal  to  t2,  as  seen  in  Figure  2.8. 

We  solve  (2.16)  for  Xi(t2) 


Xi(t2)  =Xi(0)+X2(0)t2--jNt2 


=  Xi(0)+X2(0) 


X2(0) 


-iN 


X2(0) 


=  Xi(0)  +  — x^(0). 
IV  2N  ^ 


(2.25) 


Finally  we  combine  these  solutions  to  determine  the  switching  time,  ts, 
for  any  initial  conditions  such  that  xi(0)  >  0  and  X2(0)  >  0  with  u  =  -N.  From 
Figure  2.3  we  see  that  xi(0)  for  the  ti  solution  is  the  same  point  as  Xi(t2)  from 
the  t2  solution  so  that 


and 


ts=t,+t2  = 


|x,(0) 

V  N 

|x,(tj) 

N 

x,(0)+  x^(0) 

,  2N 

1  N 

|x,(0)^l 
i  N  2 

iT) 

x,(0)  1  1 1 

'x,(0)Y 

(2.26) 


N  2K  N  J 


N 


(2.27) 


These  equations  were  derived  for  specific  initial  conditions  and  are 
only  valid  where  the  initial  conditions  lie  above  the  zero  trajectory  curves  of 
Figure  2.2.  In  order  to  expand  the  capabilities  of  ’he  control  system  for  all  initial 
conditions  we  re-derive  the  switching  time  equations  for  initial  conditions 
Xj(0)  <  0,  X2(0)  <  0,  and  u(0)  =  +N  for  the  solution 


t.  = 


X2(0) 


N 


(2.28) 


This  equation  is  valid  only  when  the  initial  values  are  below  the  zero  trajectory 
curves  of  Figure  2.2.  For  this  discussion  we  will  limit  ourselves  to  initial 
conditions  above  the  zero  trajectory  curves  and  use  the  solution  (2.27). 
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2.  Simulation  of  Analytic  Switching  Time  Controller 

Using  a  computer  simulation,  we  test  the  accuracy  of  the  solution  by 
choosing  initial  conditions  and  observing  the  response  to  our  control  effort.  We 
simulated  the  example  system  of  (2.9)  with  xi(0)  >  0,  X2(0)  >  0,  and  N  =  -1. 
The  control  effort,  u,  changes  sign  to  -N  at  the  calculated  switching  time,  t*. 
The  results  show  that  the  states  pass  through  the  origin  of  the  state  space  and 
the  error  at  the  origin  is  associated  with  the  discretization  of  the  simulation,  as 
shown  in  Figures  2.9  and  2.10.  The  switching  time  shown  in  Figure  2.10  is  the 
calculated  value.  The  control  effort  actually  switches  at  the  first  sampling  time 
following  the  caluculated  switching  time,  tj.  If  we  increase  the  sample  rate  for 
the  simulation  we  reduce  the  terminal  miss  error.  Upon  reaching  the  origin 
some  additional  control  logic  must  be  devised  to  maintain  position. 

While  the  switching  time  solution  is  a  minimum  time  solution  for  driving 
the  states  to  the  origin,  or  any  desired  values,  it  must  be  shut  off  at  tf.  The 
switching  time  solution  cannot  adapt  to  a  time  varying  situation  as  the  switching 
time  is  defined  solely  on  the  initial  conditions. 
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Figure  2.9  Simulation  of  Analytic  Switching  Time  Solution 


Figure  2.10  Control  Effort  for  Analytic  Switching  Time 
Simulation 
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C.  SOLUTION  TO  SECOND  ORDER  SWITCHING  CURVES 

An  option  for  improving  the  capabilities  of  the  second  order  controller  is  to 
remove  the  time  dependancy  of  the  control  effort.  The  control  problem  may  be 
posed  as  a  function  of  the  states  instead  of  a  function  of  time  for  some  simple 
systems.  By  defining  the  control  effort  as  only  a  function  of  the  states  the 
control  effort  is  updated  continuously  and  may  change  immediately  in  response 
to  a  change  in  the  system  parameters. 

There  are  only  two  possible  control  efforts,  -N  and  +N.  Therefore,  any 
position  in  space  will  have  either  +N  or  -N  control  effort  to  drive  it  along  its 
unique  path  to  the  origin,  and  we  can  solve  this  system  for  the  second  order 
switching  curves  that  divide  the  state  space  into  two  parts,  see  Figure  2.2. 


1.  Second  Order  Switching  Laws 
Given  the  system 


0  iTxii  ro 


^2.  .0 


(2.29) 


we  can  integrate  the  state  equations.  Setting  u  =  ±N 

X2(t)  =  I  X2(t)dt  =  J  udt  =  ±Nt  +  c,  (2.30) 

X2(t)  =  I  Xi(t)dt  =  j  X2dt  =  J±Nt  +  Cjdt  =  Nt^  +  C2t  +  C2.  (2.31) 

Choosing  the  boundary  values  so  that  xi(tf)  =  X2(tf)  =  0  we  may  rewrite  the 


equations  as 


Xi(tf)  =  ±iNtf +  X2(0)tf  +  Xi(0) 
X2(tf)  =  ±Ntf +X2(0) 


(2.32) 

(2.33) 


0  =  ±|Ntf  +X2(0)tf  +  Xj(0) 


(2.34) 
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(2.35) 

(2.36) 

(2.37) 

(2.38) 


0  =  Xj(t)  +  ^^±^^^  (2.39) 

'  N  ^  N 

0  =  x,(t)  +  i4^.  (2.40) 

N 

We  make  the  substitution 

+  X2(t)  =  -X2(t)|x2(t)l  (2.41) 

SO  that 

0  =  Xi(O-2F^2(O|x2(t)|.  (2.42) 

Equation  (2.42)  describes  the  two  parabolic  trajectories  referred  to  as 
the  zero  trajectory  curves  because  any  states  on  these  curves  will,  with  the 
appropriately  signed  control  effort,  be  driven  to  the  origin.  These  zero 
trajectory  curves  divide  the  state  space  into  two  regions  of  opposite  control 
effort.  If  the  states  are  above  the  curves  then  u  =  -N,  but  if  they  are  below  the 
curves  then  u  =  +N.  Therefore  our  control  effort,  u,  may  be  defined 

u  =  -Nsign(xi(t)-^X2(t)|x2(t)|).  (2.43) 
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With  this  switching  law  we  drive  the  states  from  any  initial  conditions  to  the 
origin  with  at  most  one  change  of  the  control  effort. 

For  example,  let  us  define  the  initial  conditions  of  the  states  to  be 
above  the  zero  trajectory  curves  so  that  Xi(0)  <  0  and  X2(0)  >  0,  as  shown  in 
Figure  2.11.  From  (2,43)  the  control  effort  is  u  =  -N,  which  drives  the  states 
along  the  parabolic  trajectory  shown  in  Figure  2.11. 

When  the  states  intersect  the  zero  trajectory  curve,  the  control  effort 
changes,  according  to  (2.43)  to  u  =  +N,  and  follows  the  zero  trajectory  curve 
into  the  origin, 

2.  Limit  Cycles 

The  control  effort,  u  =  +N  will  not  only  drive  the  states  to  the  origin,  it 
will,  in  fact,  become  confused  there.  It  is  a  point  of  indecision  and  chatter  or 
limit  cycle  motion  ensues  as  in  Figure  2.12.  The  magnitude  of  the  limit  cycle 
depends  on  the  sampling  rate  or  time  delay  for  the  control  effort.  If  the  sample 
rate  is  low,  the  states  will  penetrate  far  into  the  opposite  control  region  before 
the  control  effort  can  change,  and  the  limit  cycle  will  swing  widely  around  the 
origin.  If  the  sample  rate  is  high  then  the  states  will  not  travel  far  away  from 
the  origin  before  the  control  effort  corrects  the  direction  of  travel. 

Since  a  heavy  chatter  mode  is  not  desirable  for  many  real  systems, 
control  logic  for  reducing  or  removing  the  chatter  mode  may  be  developed, 
however,  it  will  not  be  included  in  this  paper. 
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Figure  2.11  Solution  Trajectory  From  Fixed  Initial  Condition 


Figure  2.12  Limit  Cycle  on  Second  Order  Solution  Trajectory 
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3.  Simulation  of  the  Switching  Law  Controller 

To  test  the  switching  law  (2.43),  we  simulate  the  system  of  (2.29) 
using  a  maximum  control  effort  of  N  =  1  and  initial  conditions  Xi(0)>0, 
X2(0)  >  0.  The  output  of  the  simulation,  shown  in  Figure  2.13,  demonstrates  the 
control  effort  driving  the  states  first  to  the  zero  trajectory  curve,  then  along  the 
zero  trajectory  curve  to  the  origin.  Once  at  the  origin,  the  control  effort  goes 
into  limit  cycle  as  shown  in  Figure  2.14. 
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Figure  2.13  Simulation  Using  Second  Order  Switching  Law 


Figure  2.14  Control  Effort  for  Second  Order  Switching  Law 
Solution 
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III.  APPLICATION  OF  A  SECOND  ORDER  CONTROLLER 


A.  MISSILE/TARGET  INTERCEPT  MODEL 

An  application  for  Minimum  Time  Optimal  Control  is  an  anti-missile  defense 
system,  where  the  incoming  target  has  a  speed  advantage  over  the  defensive 
missile.  In  this  case  the  missile  control  system  must  operate  in  saturation  mode. 
The  bang-bang  controller,  where  control  effort  is  either  maximum-positive  or 
maximum-negative,  has  a  faster  response  capability  than  the  standard 
Proportional  Navigation  guidance  system.  A  simple  model  of  a  missile  to  target 
engagement  can  be  described  by  Figure  3.1  where  o  is  the  line  of  sight  (LOS) 
angle  from  missile  to  target,  Ym  is  the  angle  of  the  missile  velocity,  Yt  is  the  angle 
of  the  target  velocity,  and  all  angles  are  relative  to  an  inertial  reference. 


Figure  3.1  Missile/Target  Geometry 
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It  is  assumed  that  the  target  is  on  the  final  leg  of  it's  flight  and  is  now  on  a 
straight,  non-maneuvering  trajectory.  The  control  is  to  drive  the  line  of  sight 
rate  (a)  and  it's  derivative  (a)  to  zero  in  minimum  time. 

For  our  models  we  will  be  using  only  two  dimensional  scenarios.  The 
system  dynamics  use  second  order  models  for  each  dimension.  The  target 
dynamics  are  non-maneuvering  so  that  the  acceleration,  at,  is  zero.  The  missile 
acceleration,  am,  is  assumed  perpendicular  to  the  velocity  vector  Ym.  The 
system  dynamics  are  shown  in  the  signal  flow  graph  in  Figure  3.2. 
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The  LOS  angle,  a,  and  it's  derivatives,  a  and  a,  may  be  derived  analytically 
from  the  available  states.  We  define  the  geometry  of  the  system  as  in  Figure 
3.3. 

I  JgjggJ 

I 

>% 
pZ 

II 

f 

oc 


Figure  3.3  Geometry  for  Angle  Definitions 
with 

Rt  =  Range  of  the  Target. 

Rm  =  Range  of  the  Missile. 

Vt  =  Velocity  of  the  Target. 

Vm  =  Velocity  of  the  Missile, 
at  =  Acceleration  of  the  Target, 
am  =  Acceleration  of  the  Missile. 
Vc  =  Closing  Velocity. 

so  that 


CT  =  tan 


f  R  ^ 

*'tmy 


(3.1) 
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o  = 


R  .v  — R  V 

^tmx  ^tmy  ^tmx 


(3.2) 


and 


^tmx  ‘  ^tmy  ^^tmx  '  ^tmy  ^tmy  '^tmx  ,  ^^c(^tn«  ^tmy  ^tmy  ^tmx) 


^tmy  *^tmy  ’•imx  ,  u»»a  uuj  uhj  mhw  ^ 

a  =  - -^5 - + - ^3  (3.3) 

where 

V,=-R  (3.4) 

In  an  actual  application  Kalman  or  Luenberger  Observers  may  be  used  to 
obtain  estimates  of  a  and  c. 


B.  APPLICATION  OF  THE  SWITCHING  LAW 

The  switching  law  from  (2.43)  is  applied  to  our  scenario  where  6  and  d 
form  the  state  space  of  this  simulation,  and  is  implemented  as 


u  =  N  •  sign 


a  +  - 


G_o 

2N 


(3.5) 


The  sign  convention  of  (2.43)  was  changed  to  conform  to  the  geometry  of  our 
scenario. 

We  set  the  simulation  with  the  initial  conditions 


Rmx(O)  =  0  ft 

Vn,x(0)  =  2(XX)  ft/sec 
amx(O)  =  0  ft/sec2 
R«(0)  =  10000  ft 
V„(0)  =  2000  ft/sec 
atx(O)  =  0  ft/sec2 
tfinai  =  2.25  sec 


Rn,y(0)  =  0  ft 
v„y(0)  =  0  ft/sec 
amy(O)  =  0  ft/sec2 
Rty(O)  =  1000  ft 
V,y(0)  =  0  ft/sec 
aiy(O)  =  0  ft/sec2 
dt  =  0.01  sec 
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Running  the  simulation  we  find  that  the  controller  does  drive  6  and  a  to 
zero,  and  maintains  them  about  the  origin  until  intercept  is  reached  (Figure  3.4). 
Once  the  states  reach  the  origin,  the  system  controller  goes  into  chatter  mode, 
see  Figure  3.5,  and  the  trajectory  chatters  back  and  forth  about  the  desired 
intercept  path. 
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Figure  3.4  Missile/Target  Intercept  Simulation  Using  Second 
Order  Controller 
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Figure  3.5  Control  Effort  for  Second  Order  Missile/Target 
Simulation 
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C.  COMPARISON  WITH  PROPORTIONAL  NAVIGATION 
CONTROLLER 

The  system  dynamics  for  the  proportional  navigation  controller  simulation  is 
basically  the  same  as  that  in  the  previous  section,  except  that  the  control  effort , 
u,  is  proportional  to  d,  which  is  obtained  with  an  estimator,  shown  in  Figure  3.6, 

•  A 

where  |3  =  d.  We  limit  the  acceleration  so  that  u  is  bounded  by  ±N.  With  only 
the  change  in  the  controller  we  ran  the  simulation  for  the  same  initial  conditions 
with  the  results  shown  in  Figures  3.7  and  3.8. 

This  Proportional  Navigation  system  is  effective  at  intercepting  the  target 
only  when  the  missile/target  geometry  and  kinematics  are  sufficient  that  the 
missile  has  time  to  maneuver.  The  delay  in  saturation  of  the  control  effort 
caused  the  missile  to  turn  too  slowly  so  that  the  Proportional  Navigation 
Controlled  system  was  unable  to  maneuver  quickly  enough  to  intercept  the 
faster,  although  non-maneuvering,  target. 
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Figure  3.6  System  Dynamics  of  Proportional  Navigation 
Missile/Target  Model 
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Figure  3.7  Proportional  Navigation  Missile/Target  Simulation 


Figure  3.8  Control  Effort  for  Proportional  Navigation 
Missile/Target  Simulation 
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IV.  THIRD  ORDER  CONTROLLER 


A.  FORWARD  TIME  SYSTEM 

The  second  order  controller  is  effective  at  intercepting  a  target  that  has 
speed  advantage.  Previously  we  have  controlled  the  system  by  driving  c  and 
d  to  zero,  and  thereby  maintaining  a  constant  LOS  angle,  o,  until  impact.  But 
what  about  controlling  the  LOS  angle  itself?  There  are  situations  in  which  we 
would  like  to  be  able  to  attack  a  target  from  a  particular  angle,  or  perhaps  have 
multiple  missiles,  each  with  it's  own  pre-defined  attack  angle. 

Consider  a  point  defence  system  in  which  the  missile,  launched  from  some 
point  away  from  the  target's  final  trajectory,  would  try  to  position  itself  in 
minimum  time  onto  a  "head-on"  collision  course  with  the  target,  as  shown  in 
Figure  4.1.  In  such  a  situation  we  would  use  a  third  order  minimum  time 
controller  to  drive  a,  a,  and  a  to  zero. 

1.  System  Definition 

Our  third  order  example  is 


"0 

1 

0' 

'O' 

X  = 

0 

0 

1 

x-t- 

0 

0 

0 

0 

1 

Minimizing  the  Hamiltonian  and  solving  the  system  we  find 

H  =  1  +  PiX2+P2X3  +  P3U 


and 


u  =  -Nsign(p3) 


(4.1) 


(4.2) 

(4.3) 
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Figure  4.1  Application  of  Third  Order  Minimum  Time  Controller 
where 


so  that 


'0  0  O' 

0  -1  0  p 

0  0  -1 


Pi  =Cj 


(4.4) 


(4.5) 


P2  =  -Cit  +  C2 

(4.6) 

P3  =^Cit^-C2t  +  C3. 

(4.7) 
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From  the  adjoint  solution  the  control  may  switch  sign  no  more  than 
twice.  Again  tracing  this  problem  in  negative  time  we  may  follow  the  zero 
trajectory  curves  out  from  the  origin  with  control  efforts  of  ±N.  These  curves 
are  now  in  three  dimensional  space  residing  in  their  own  plane  which  intersects 
the  origin.  Intersecting  these  zero  trajectory  curves  are  an  infinite  number  of 
curves  making  a  surface,  and  leading  off  from  this  surface  the  infinite  number  of 
trajectories  lead  to  the  initial  conditions.  Therefore  in  forward  time  we  may 
start  with  an  initial  condition  such  that  u  =  +N  will  drive  the  system  to  intersect 
with  the  surface  at  tjwi  as  shown  in  Figure  4.2.  Switching  the  control  effort  to 
u  =  -N  will  drive  the  system  along  the  surface  to  intersect  with  the  zero 
trajectory  curve  at  tsw2  where  the  control  effort  switches  again.  Finally  u  =  +N 
will  drive  the  system  to  the  origin. 


Figure  4.2  Three  Dimensional  Switching  Curves 
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2.  Third  Order  Switching  Curves 

We  have  determined  from  the  second  order  system  that  solving  for  the 
control  switching  times  is  not  as  widely  applicable  as  defining  the  control  as  a 
function  of  the  states;  so  we  now  move  on  to  defining  the  third  order  switching 
curves. 

Starting  with  the  state  equations  (4.1)  we  discretize  and  expand  the 
equations  so  that 


(j)  =  e"^'  = 


1  + At  + 


A^t^  A^t^ 


2! 


3! 


‘+... 


’1 
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o' 

'0 
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t^+. 


1  t2 


1  t  ^t 


0  1 
0  0 


A3  is  the  zero  matrix  so  it  and  all  higher  terms  of  A  are  dropped. 


(4.8) 


A=  re^^‘~^^Bdx=  f 
Jo  Jo 


1  t-T  ^(t-T)^ 
0  1  t-T 

0  0  1 


'o' 

0 

1 

dx 


■it^x-itx^  +  Jx^' 

t 

6  '■ 

t-x 

dx  = 

tx--|x^ 

= 

J-t^ 

2  *■ 

1 

X 

0 

t 

(4.9) 


Combining  these  we  obtain 
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(4.10) 


'i  t  itn  rit^' 

x(t)=  0  1  t  x(0)+  u(0) 

0  0  1  t 

or  in  scalar  equations 

Xi(t)  =  Xi(0)  +  tX2(0)+it^X3(0)+^t^u(0)  (4.11) 

X2(t)=  X2(0)  +  tX3(0)  +  -it^u(0)  (4.12) 

X3(t)  =  X3(0)  +  tu(0).  (4.13) 

B.  NEGATIVE  TIME  SYSTEM 

The  third  order  system,  being  piecewise  continuous,  is  easily  broken  into 
several  simple  boundary  value  problems.  In  order  to  solve  the  systems  of 
curves  we  must  determine  some  boundary  values  for  the  intersections  of  these 
families  of  curves.  We  start  at  the  origin  and  run  the  system  in  negative  time  to 
determine  our  other  boundary  conditions. 

In  negative  time  we  have 

'0-10  1  [O' 

x=  0  0  -1  x+  0  u.  (4.14) 

0  0  oj  [-1 

Discretizing  we  find 

fl  -t  it^' 

<(,  =  e^‘=  0  1  -t  (4.15) 

0  0  1 

and 
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(4.16) 


A=  = 

Jo 


6  '■ 

it^ 
2  ^ 

-t 


SO  that 


"l  -t  It^' 

6  ^ 

X(t)  = 

0  1 

-t 

x(0)+ 

it^ 

2  ‘ 

0  0 

1 

-t 

u(0), 


(4.17) 


or  in  scalar  equations 


Xi(t)  =  Xi(0)-tX2(0)+it^X3(0)-it^u(0) 

(4.18) 

X2(t)=X2(0)-tX3(0)  +  it^u(0) 

(4.19) 

X3(t)  =  X3(0)-tu(0). 

(4.20) 

1.  Solving  for  Negative  Time  Boundary  Points 

To  develop  a  complete  solution  we  solve  the  equations  for  several 
different  points  along  the  zero  trajectory  curve.  Setting  u  =  -N  and  traveling  out 
along  the  zero  trajectory  curve  for  1  second  we  find 


x,(l)  =  -iut’  =  -i{-N)  =  iN 

(4.21) 

X2(l)  =  iut"=^(-N)  =  -iN 

(4.22) 

X3(l)  =  -ut  =  -(-N)  =  N. 

(4.23) 

Similarly  we  run  the  system  in  negative  time  for  2  and  3  seconds  to 
obtain  other  boundary  points  as  shown  in  Figure  4.3,  and  listed  in  TABLE  1. 
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Figure  4.3  Boundary  Values  in  Negative  Time  Solution 


NEGATIVE  TIME  BOUNDARY  CONDITIONS,  u  =  -N 


u  =  -N 

t  =  2 

t  =  3 

xi(t) 

jN 

fN 

X2(t) 

-2N 

-fN 

X3(t) 

2N 

3N 

TABLE  1. 


The  control  effort  may  also  be  u  =  +N,  traveling  away  from  the  origin 
on  the  other  zero  trajectory  curve  giving  the  boundary  conditions  in  TABLE  2. 
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NEGATIVE  TIME  BOUNDARY  CONDITIONS,  u  =  +N 


u  =  +N 

t=  1 

xi(t) 

-iN 

X2(t) 

In 

X3(t) 

-N 

-2N 


2N 


TABLE  2. 

We  may  now  solve  the  equations  for  the  family  of  curves  that  intersect 
the  zero  trajectory  curves.  Specifically,  we  solve  for  the  equation  of  the  one 
curve  that  travels  from  some  point  xi(0),  X2(0),  X3(0),  to  the  point  Xi(t)  =  -^N, 


X2(t)  =  ^N,  X3(t)  =  N.  Since  the  control  effort  on  the  zero  trajectory  curve  for 
this  intersection  point  is  u  =  -N,  the  control  effort  for  the  curve  we  are  solving 
for  must  be  u  =  +N,  and  the  forward  time  equations  may  be  written  as 


Xi(t)  =  iN  =  Xi(0)+tX2(0)  +  it%(0)  +  it^N 
X2(t)  =  -iN=X2(0)+tX3{0)  +  it^N 
X3(t)=N  =  X3(0)+tN 
Solving  (4.26)  for  t  we  find 

N 


(4.24) 

(4.25) 

(4.26) 


(4.27) 


Substituting  (4.27)  into  (4.24)  and  (4.25),  and  after  a  bit  of  algebra  we  obtain 

KT  _  ^3  (A 


-N  =  x2(0)- 


0  =  Xi(0)+  X2(0)- 


X2(0)X3(0)  X32(0)  X33(0) 


(4.28) 

(4.29) 


'  ^  N  2N  3N2 

Using  the  other  boundary  values  from  TABLE  1  and  TABLE  2  as 
well,  we  may  generate  a  family  of  equations  representing  both  sides  of  the  zero 


trajectory  curves,  and  different  points  of  intersection  along  the  curves  (see 
TABLE  3).  The  control  effort,  u,  in  TABLE  3.  is  the  control  effort  for  the  zero 
trajectory  curve  that  is  intercepted.  The  time,  t,  is  the  time  out  from  the  origin 
to  the  intercept  point  for  the  negative  time  system. 

To  combine  all  the  equations  from  TABLE  3  into  one  solution  we 

define 


RvMILY  OF  SOLUTION  CURVES 


u  =  -N 

-N  =  x2(0)-^^^^ 

2N 

t  =  1 

0-x(0)lx  (0)  ’‘2(0)>‘3(0)  ’‘3^(0)  Xs’CO) 
0-x,(0)  +  x,(0)  ^ 

u  =  -N 

-4N  =  x2(0)--^^^-^ 

2N 

t  =  2 

0  =  x,(0)  +  2x2(O)  M(0)X5(0)_x/(0)^x,’(0) 

^  ^  N  N  3N^ 

u  =  -N 

-9N  =  x2(Q)--^^^-^ 

2N 

t  =  3 

*  ^  ^  N  2N  3N^ 

u  =  +N 

N  =  x2(0)+^^2!i22 

2N 

t=  1 

^  ^  N  2N  3N^ 

u  =  +N 

4N  =  x2(0)+^^ 

2N 

t  =  2 

0  =  x,(0)3-2x2(0)+ 

‘  ^  ^  N  N  3N^ 

u  =  +N 

9N  =  x2(0)-t-^^-^ 

2N 

t  =  3 

^  ^  N  2N  3N2 

TABLE  3. 


C.  THIRD  ORDER  SWITCHING  LAW 


Using  (4.30),  (4.31)  and  (4.32)  we  can  break  up  space  into  areas  of 
opposite  control  effort.  The  control  effort  is  defined  by  which  of  these  volumes 
contain  the  states.  Therefore,  the  third  order  switching  law  is  defined  as 


u  =  -N  sign  Xi(0)  + 


^3^(0)  ,  X2(0)  X3(0) 


(4.33) 


D.  THIRD  ORDER  CONTROLLER  SIMULATION 

A  simulation  of  this  third  order  model  shows  that  the  third  order  switching 
law,  (4.33),  drives  the  states  to  the  origin  with  only  2  changes  is  the  control 
effort.  In  application  some  means  of  shutting  off  the  control  effort  will  be 
required  for  the  system  to  avoid  entering  chatter  mode  upon  reaching  the  origin 
(see  Figure  4.4). 

An  examination  of  the  functions  (4.30),  (4.31),  and  (4.32)  in  Figure  4.5 
shows  that  f  is  a  smoothly  increasing  curve  until  the  states  intercept  the 
switching  surface.  Once  the  states  are  following  the  surface,  f  remains  a 
constant  value.  When  the  states  reach  the  zero  trajectory  curve,  f  becomes 
approximately  zero. 


40 


Control  Functions 


Figure  4.4  Simulation  of  a  Third  Order  Minimum  Time 
Controller  From  a  Fixed  Initial  Condition 


Figure  4.5  Control  Laws  For  Third  Order  Solution 
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V.  MISSILE  ADJOINTS 


Given  a  time  varying  system 

x  =  A(t)x  +  B(t)u  (5.1) 

y  =  C(t)x.  (5.2) 

it  can  be  shown  that  the  impulse  response  of  the  adjoint  system 

p  =  A(t,-tf  p  +  C(t,-tfr  (5.3) 

y  =  B(l(-t)’^p  (5.4) 

is  the  response  of  the  original  system,  at  time  t,  to  an  impulse  applied  at  time 
tf  - 1  before  t.  [2]  For  a  time  invariant  system  the  transfer  function  for  the 
original  system  is  identical  to  the  transfer  function  of  the  adjoint  system,  i.e., 
they  are  self-adjoint. 

A.  CONSTRUCTION  OF  AN  ADJOINT 

Given  the  mathematics  of  the  adjoint  method,  we  now  need  to  define  some 
rules  to  create  and  understand  the  adjoint  system  with  real  systems  and  block 
diagrams.  [3] 

1.  Rule  1:  Convert  All  System  Inputs  to  Impulses 

In  constructing  the  adjoint  it  is  necessary  that  all  system  inputs  be 
impulsive.  Since  this  may  not  be  the  case  with  the  block  diagram,  all  system 
inputs  and  initial  conditions  must  be  converted  to  impulsive  inputs  via  block 
manipulations  and  extensive  use  of  integrators..  Figure  5.1  shows  that  step 
inputs  and  initial  conditions  are  equivalent  to  the  output  of  impulse  driven 
integrators. 
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2.  Rule  2:  Replace  t  With  tf  -  t  in  the  Arguments  of  All  Time 
Varying  Coefficients 

In  many  linear  systems  it  is  possible  to  express  a  gain  as  a  function  of 
time.  The  adjoint  system  operates  in  negative  time  and  Figure  5.2  shows  the 
conversion  of  functions  of  time  into  the  adjoint  domain. 


Orieinal  Sv.stem 

Adioint  Svstem 

Time 

Varying 

K(t)  =  at+b 

K(l>=  '  ' 

a(t-  tf )  +  b 

K(t-tf)=a(t-tf)  +  b 

Gain 

K(t-tf)=— 

at  +b 

Figure  5.2  Convert  Functions  of  Time  to  the  Adjoint  Domain 
3.  Rule  3:  Reverse  the  Direction  of  all  Signal  Flow 


The  direction  of  all  signal  flow  must  be  reversed,  redefining  nodes  as 
summing  junctions  and  visa  versa.  Notice  that  all  system  outputs  become 
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inputs  and  all  system  inputs  become  outputs.  This  last  rule  allows  for  the  simple 
graphical  creation  of  the  adjoint  system  by  first  drawing  the  block,  or  flow, 
diagram  of  the  original  system,  then  redrawing  it  with  the  arrows  reversed. 
Figure  5.3  shows  some  examples  f  converting  nodes  to  summing  junctions  and 
visa  versa. 


Original  System 


Adjoint  System 


G(s) 

— 


m 

— >► 


<  •- 


G(s) 


Yb(t-tf) 
— < - 


Figure  5.3  Converting  Nodes  to  Summing  Junctions 


B.  DEVELOPMENT  OF  A  SECOND  ORDER  ADJOINT  SYSTEM 
This  adjoint  formulation  lends  itself  well  to  analyzing  some  optimization 
problems,  those  where  tf  is  free  and  the  terminal  state  is  constrained  to  a  point, 
curve  or  surface.  From  Pontryagin  the  control  is  a  function  of  the  homogeneous 
adjoint  (see  (4.2)  and  (4.3)).  The  missile  adjoint  form  gives  also  a  particular 
solution  of  the  adjoint  (i.e.  system  impulse  response). 

The  adjoint  allows  one  to  generate  optimum  solutions  (i.e.  switching 
surfaces)  for  all  possible  initial  conditions.  The  forcing  impulses  passed  through 
an  integrator  gives  our  saturation  type  of  optimal  control  (±N). 

As  an  example  for  the  application  of  the  method  of  adjoints  we  will  apply 
our  adjoint  rules  to  our  second  order  controller.  The  system  described  by  (2.9) 
has  a  time  dependent  control  input  where 
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(5.5) 


u(t)  = 


t<ts 

t>ts 


and,  from  our  second  order  example. 


ts  = 


x(Q)  ^  l|<x(0)Y  ^  x(0) 
N  2I  N  J  N  ' 


(5.6) 


1.  Apply  Rule  1  to  the  Second  Order  Example 

For  this  example  we  will  define  N  =  1 .  We  first  draw  out  the  original 
signal  flow  diagram,  showing  all  the  inputs  and  initial  conditions.  In  order  to 
apply  the  impulsive  inputs  we  expand  the  second  order  system  to  the  third 
order  state  equations 


‘0 

1 

o' 

'O' 

X  = 

0 

0 

1 

x  + 

0 

_0 

0 

0 

1 

(5.7) 


y  =  [l  0  0]x.  (5.8) 

Converting  all  the  inputs  to  impulsive  inputs  and  initial  conditions  we  get  the 
flow  diagram  in  Figure  5.4. 
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Figure  5.4  Application  of  Rule  1  on  Second  Order  Example 
2.  Apply  Rule  2  to  the  Second  Order  Example 

We  replace  t  with  t  -  tf  in  the  arguments  of  all  time  varying  coefficients 
to  obtain  the  flow  diagram  in  Figure  5.5. 
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3.  Apply  Rule  3  to  the  Second  Order  Example 

Finally  we  reverse  the  direction  of  the  signal  flow,  redefining  nodes  as 
summing  junctions  and  visa  versa,  thereby  changing  the  system  inputs  to 
outputs  and  the  system  outputs  to  inputs  as  shown  in  Figure  5.6. 


The  mathematical  definition  of  the  Adjoint  System  is 

P  =  a’^P+C'^u  (5.9) 

y  =  B^p.  (5.10) 

Using  the  A,  B,  and  C  matrix  from  (5.7)  and  (5.8)  we  obtain 


'pl' 

'0 

0 

o' 

Pl" 

T 

P2 

= 

1 

0 

0 

P2 

+ 

0 

.P3. 

_0 

1 

0 

-P3. 

0 

y  =  [0  0  l]p 

which  corresponds  to  Figure  5.6. 

C.  SIMULATION  OF  THE  SECOND  ORDER  ADJOINT 
1.  Forward  Time  Second  Order  Simulation 


(5.11) 

(5.12) 


In  the  forward  time  second  order  simulation  we  set  N  =  1,  and  the 
initial  conditions  x(0)  =  2  and  x(0)  =  0.  Using  impulses  through  integrators  we 


apply  the  initial  control  effort  at  t  =  0.  A  second  impulse,  opposite  in  sign  and 
twice  in  magnitude,  is  applied  at  the  switching  time,  as  defined  in  (5.6),  driving 
the  states  through  the  origin. 

The  simulation  length  is  4  seconds  with  a  sample  step  of  .05  seconds. 
The  results  are  shown  in  Figure  5.7.  Starting  at  the  initial  conditions  the  system 
is  driven  through  the  origin  with  a  minimum  miss  distance  of  0.003344  at  a  time 
of  2.83  sec. 

2.  Adjoint  Solution  for  Second  Order  System 

In  the  adjoint  domain  the  system  will  travel  in  negative  time  starting  at 
the  origin  and  traveling  outward  to  the  initial  conditions  for  the  forward  time 
system.  Because  this  is  a  time  invariant  system  the  trajectory  for  the  adjoint 
solution  will  be  the  same  as  the  forward  time  system  but  in  the  opposite 
direction.  From  our  second  order  example,  the  optimum  switching  in  negative 
time  from  the  terminal  state  at  the  origin  is 


If  ““  1-  . 

f  *  N 


The  trajectory  constraint  yields 


(5.13) 


X2(ts)|x2(ts)| 

2N 


(5.14) 


A  reverse  impulse  of  twice  the  magnitude  is  applied  and  the  trajectory 
proceeds  out  to  the  desired  initial  conditions  given.  The  negative  time  impulse 
response  is  thus  prescribed  by  using  the  switching  times 


tf-ts  =  ^VNx,(0)  +  ix2(0)2 


(5.15) 
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tsw  -to  =-^VN  x,(0)+ixj(0)^ 


(5.16) 


The  parameters  for  the  adjoint  solution  are  the  same  as  for  the 
forward  time,  but  the  initial  values  for  the  states  are  at  the  origin.  The  output 
of  the  adjoint  system,  according  to  (5.12),  is  ps,  so  the  phase  plot  of  Figure  (5.8) 
plots  -p2  and  p3.  The  adjoint  solution  traces  almost  the  same  path  as  the 
forward  time  solution,  switching  at  tj  =  1.582  sec.,  and  missing  the  point  (2,1)  by 
0.002733  at  tf  =  2.83  seconds.  By  changing  the  initial  sign  of  the  control  effort, 
and  adjusting  the  switching  time,  we  could  drive  the  adjoint  system  to  any 
desired  point  in  the  phase  plane,  and  therefore  know  the  optimal  solution  for  the 
forward  time  system. 


49 


Figure  5.8  Simulation  of  the  Adjoint  Solution 
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D.  MISSILE  SIMULATION  WITH  ADJOINTS 


The  method  of  adjoints  can  be  extremely  useful  when  dealing  with  time 
varying  systems.  By  using  the  adjoint  solution  we  are  able  to  see  how  the 
forward  time  system  behaves  for  all  times. 

1.  Missile/Target  Model 

We  will  simplify  our  scenario  of  missile/target  engagement  in  order  to 
present  an  uncluttered  example  of  the  method.  The  target  will  have  zero  x- 
velocity  and  constant  y-velocity,  and  will  start  on  the  x-axis  at  a  distance,  R, 
from  the  origin  (see  Figure  5.9).  The  missile  will  start  at  the  origin  with  a 
constant  x-velocity  by  zero  initial  y-velocity. 


Figure  5.9  Geometry  of  Missile/'iarget  Adjoint  Solution 

Since  R  is  targe  and  yt  -  ym  is  small  we  use  the  small  angle 
approximation  and  define 


a  =  tan' 


yt-ym 


R 


R 


The  closing  velocity  is  constant  so  we  note  that 

R  =  Vc(tf-t) 


(5.17) 

(5.18) 
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where  tf  is  the  length  of  the  simulation. 

2.  Signal  Flow  Diagram 

Because  the  x-velocities  are  constant  we  need  only  to  model  the  y- 
dimension.  We  develop  the  flow  diagram  in  Figure  5.10  from  the  geometry  of 
Figure  5.9.  We  use  proportional  navigation  control  and  impulsive  inputs  are 
used  for  initial  conditions. 


The  estimator  used  was  developed  from  the  mechanics  and  dynamics 
of  a  seeker  head  system,  however,  it  can  be  shown  to  be  equivalent  to  a 
Kalman  Estimator  or  Luenberger  Observer.  This  estimator  in  Figure  5.11  has  a 
time  constant  of  0.1  seconds.  The  primary  interest  in  this  model  is  to  determine 
the  miss  distance  at  the  final  time,  t  =  tf,  so  the  output  of  the  system  is  ymiss- 
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3.  Forward  Time  Simulation 
We  define  our  states  to  be 


'ym' 

X2 

ym 

X3 

X4 

p 

X5 

yt 

-X6. 

.yi  . 

and  the  system  state  equations  are  then 


X  = 


0 

0 

0 

-100 

tf-t 

0 

0 


1 

0 

0 

0 

0 

0 


0 


0 


0  nV. 


0 


1 


-100  -20 

0  0 


0 


0 


0 

0 

0 

100 

tf-t 

0 

0 


0 

0 

0 

0 

1 

0 


x-f 


u 


y  =  [-l  0  0  0  1  0]- 


(5.19) 


(5.20) 


(5.21) 
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Since  the  x-velocities  are  constant  we  know  approximately  when  the 
point  of  closest  approach  will  occur,  and  therefore  about  how  long  to  run  the 
simulation.  The  LOS  angle,  a,  is  a  function  of  t  and  tf,  so  that  different  values  of 
tf  will  result  in  changes  to  the  control  effort  and  system  response.  Because  this 
state  matrix.  A,  is  time  dependent ,  it  must  be  redefined  at  each  time  step  of  the 
simulation.  We  define  our  proportionality  constant,  n,  to  4,  the  closing  velocity, 
Vc,  to  5000  ft/sec.,  so  that  with  tf  =  4  sec.  the  initial  range  is  20,000  ft.  The 
output  of  the  simulation  is  shown  in  Figure  5.12. 

To  study  the  effect  of  different  values  of  tf,  or  to  locate  the  optimal 
value,  we  may  have  to  run  the  forward  time  simulation  many  times,  which 
could  be  a  very  tedious  process,  especially  since  we  are  only  interested  in  the 
final  value  of  t  =  tf. 

4.  Adjoint  Solution 

An  alternative  to  running  the  forward  time  system  over  and  over  again 
would  be  to  run  the  adjoint  solution  once  to  generate  the  final  values  of  the 
family  of  forward  time  solutions.  This  time  we  generate  the  adjoint  system 
using  matrix  algebra  instead  of  the  graphical  method.  The  solutions  to  the 
adjoint  system  are,  from  (5.9)  and  (5.10) 
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0  0 
1  0 


0 

0 


-100 

tf-t 

0 


0 

0 


0 

nY; 


0  -100  0 

1  -20  0 


0 

0 


0 

0 


0 

0 


tf-t 
0  1 


0 

0 

0 

0 


-1 

0 

0 

0 

1 

0 


•u 


(5.22) 


y  =  [0  0  0  0  0  l]  p.  (5.23) 

This  system  generates  a  curve  whose  value  at  any  time,  ta,  is  the  final 
value  of  the  forward  time  system  where  tf  =  ta.  Figure  5.13  shows  four  curves 
from  the  forward  time  solution  corresponding  to  tf  =  1,2,3,  and  4  sec.  Each 
curve  end  exactly  on  the  adjoint  solution  curve  at  the  corresponding  adjoint 
time. 
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Figure  5.13  Adjoint  Solution  Curve 
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VI.  CONCLUSIONS  AND  COMMENTS 


We  have  developed  the  second  order  minimum  time  optimal  controller  and 
applied  it  to  the  fast  reaction  missile  defense  problem  with  both  the  closed  form 
analytic  solution,  and  the  open  form  switching  curves  solution.  In  comparing 
our  open  form  solution  with  a  standard  proportional  navigation  controller 
solution  we  demonstrate  an  increased  maneuverability  enabling  our  missiles  to 
defend  against  faster,  more  capable  threats.  The  accuracy  at  intercept  is  a 
function  of  the  control  logic  used  to  shut  off  the  control  effort  when  desired 
conditions  are  met.  This  is  a  subject  that  should  be  explored  in  future  projects. 

We  introduced  a  third  order  minimum  time  controller  which  promises  to  not 
only  improve  reaction  time  and  maneuverability,  but  also  presents  us  with  the 
ability  to  control  and  define  a  desired  attack  angle  for  an  improved  destructive 
potential.  A  model  for  a  practical  third  order  controller  should  be  developed 
and  evaluated. 

Having  introduced  the  method  of  adjoints  and  shown  some  of  its  functions, 
we  hope  to  stimulate  some  more  practical  applications  of  this  technique.  We 
are  excited  at  the  possibilities  of  using  the  method  of  adjoints  in  determining 
optimal,  closed  form  solutions  to  forward  time  problems  at  speeds  fast  enough 
for  practical  implementation. 
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APPENDIX  A.  PROGRAM  CODE 


All  of  the  simulations  for  this  project  were  run  on  both  IBM-ATi  class  and 
Macintosh  IP  computers  using  the  matrix  manipulation  language  MATLAB^. 
For  IBM- AT  compatibles  MATLAB,  version  3.5f  was  used,  and  on  the 
Macintosh  II  computers  version  MAC  II-MATLAB,  version  1.1b.  This 
appendix  contains  the  source  code  for  all  of  the  simulations  and  functions 
written  in  support  of  this  project. 

Only  a  limited  background  in  programming  is  required  for  understanding 
these  files.  While  MATLAB  is  similar  to  FORTRAN,  MATLAB's  control 
structures  are  much  less  complex,  and  with  matrix  manipulation  built  into  the 
system,  vector  definition  and  storage  are  greatly  simplified.  Comments  are 
started  by  the  percent  sign  (%)  and  continue  to  the  end  of  the  line.  Ellipsis  (...) 
at  the  end  of  a  line  indicate  the  continuance  of  the  logical  line  onto  the  next  line 
of  code. 

Each  file  will  begin  on  a  new  page  to  assist  those  who  are  interested  in 
examining  or  reproducing  the  code.  Although  an  analysis  of  these  files  is  not 
necessary  to  understand  this  report,  readers  are  encouraged  to  examine  them 
closely  for  further  information. 

The  code  is  presented  in  the  order  of  usage  in  the  main  text  with  all 
supporting  functions  grouped  with  the  main  program  of  interest. 


'  IBM  and  IBM-AT  are  registered  trademarks  of  IBM. 

^  MAC  II  is  a  registered  trademark  of  APPLE. 

3  MATLAB  is  a  registered  trademark  of  The  MathWorks,  Inc.  [4] 
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1.  BB2NDST.M 


%  BB2NDST.M  11  Mar.  1991 

%  BB2NDST.M  is  a  simulatior.  of  the  2nd  order  bang-bang  controller  with  an  analytic 
%  solution  for  the  Switching  T.me  of  the  control  effort.  The  signs  of  the  control  effort 
%  must  be  matched  with  the  initial  conditions  to  have  convergence. 

%  written  by  Colin  R.  Cooper 

%  Define  the  State  Equations. 

A  =  [0  1;0  0]; 

B  =  [0  1]'; 

C  =  [10]; 
xl0  =  2; 
x20  =  1 ; 

N  =  1 ;  %  Maximum  control  effort 
aN  =  abs(N); 

Tf  =  4.17; 
dt  =  .005; 

[Phi, Dei]  =  c2d(A,B,dt); 

%  Create  the  storage  vectors, 
kmax  =  Tf/dt  +  1; 

X  =  zeros(2,kmax); 
y  =  zeros(l,kmax); 
time  =  zeros(l,kmax); 

x(:,l)  =  [xl0;x20];  %  Set  initial  conditions  for  x. 

%  Define  the  Switching  Time  for  the  control  effort, 
tsf  =  abs(x20)/aN  +  sqrt(abs(xlO  +  x20*abs(x20)/(2*aN))/aN); 

%  Begin  simulation  loop. 


%  xl  initial  condition. 

%  x2  initial  condition. 

%  Length  of  simulation. 

%  Time  increment  for  simulation. 
%  Discretized  System. 
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for  i  =  1  rkmax  -  1 
if  tiine(i)<tsf 
u  =  -N; 
else 

u  =  N; 
end 

x(:,i+l)  =  Phi*x(:,i)  +  Del*u; 
y(l,i+l)  =  C*x(:,i+1); 
time(i+l)  =  time(i)-Kit; 
end 

%  Plot  the  output  of  the  simulation. 
clg,plot(x(l,;),x(2,:),'-w’),grid 
title('Phase  plot  -  Switching  Time') 
xlabel(’Xr),ylabel('X2') 


2.  BB2NDSL.M 


%  BB2NDSL.M  11  Mar.  1991 

%  is  a  simulation  of  the  2nd  order  bang-bang  controller  using  a  Switching  Law  for  the 
%  control  effort. 

%  written  by  Colin  R.  Cooper 

A  =  [01;00];% 

B  =  [0  U'; 

C  =  [10]; 
xlO  =  2; 
x20=  1; 

N  =  1 ;  %  Maximum  control  effort 
Tf  =  4.5; 
dt  =  .01; 

[Phi, Del]  =  c2d(A,B,dt); 

%  Create  storage  vectors, 
kmax  =  Tf/dt  +  1 ; 

X  =  zeros(2,kmax); 
y  =  zeros(l,kmax); 
u  =  zeros(l,kmax); 

x(:,l)  =  [xl0;x20];  %  Set  initial  conditions  for  x. 

%  Begin  simulation  loop, 
for  i  =  1  :kmax  -  1 

u(i)  =  -N*sign(x(l,i)  +  .5*x(2,i)*abs(x(2,i))/N); 
x(:,i+l)  =  Phi*x(:,i)  +  Del*u(i); 

(l,i+l)  =  C*x(:,i+l); 

end 


Define  the  State  Equations. 

%xl  initial  condition. 

%  x2  initial  condition. 

%  Length  of  simulation. 

%  Time  increment  for  simulation. 
%  Discretize  the  System. 
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%  Plot  the  output  of  the  simulation. 
clg,plot(x(l,:),x(2,:)),grid 
title('Phase  plot  -  Switching  Law') 
xlabel(’Xl'),ylabel('X2') 


3.  SIM2SL.M 

%  SIM2SL.M  2nd  Order  Switching  Laws  Control  08  Apr.  1991 

%  Simulation  of  the  missile/target  simulation. 

% 

%  ->  Uses  analytic  values  for  sigma-dot, sigma-ddot  for  calcultion  of  the  control  effort 
%  using  the  2nd  order  switching  curves. 

%  ->  The  Target  makes  no  evasive  maneuvers. 

%  ->  Calculates  the  Quantization  Error  based  on  the  average  velocities  for  the  crossover 
%  endpoints. 

%  ->  Calls  INTERP.M  function  which  takes  the  states  at  the  crossover  endpoints,  creates 
%  100  point  straight  lines  to  connect  the  points,  and  determine  the  minimum  miss 

%  distance. 

%  ->  Allows  delay  time  for  missile  control. 

%  ->  Target  is  now  on  level  flight  with  Beta  =  0. 

%  written  by  Colin  R.  Cooper 


%  Define  states. 

N  =  1000; 

Am  =  [0  1  0  0;0  0  0  0;0  0  0  1;0  0  0  0]; 
Bm  =  [0  0;1  0;0  0;0  1]; 

At  =  [0  1  00;0000;000  1;00001; 
Bt  =  [0  0;10;0  0;01]; 

Tf=2.25; 
dt  =  .01; 

%  Descretize  the  states. 

[Phim,Delm]  =  c2d(Am,Bm,dt); 
[Phit,Delt]  =  c2d(At,Bt,dt); 


%  Maximum  control  effort. 
%  Missile  State  Equations. 

%  Target  State  Equations. 

%  Total  time  of  simulation. 
%  Sample  step  size. 

%  Discrete  Missile  System 
%  Discrete  Target  System 
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%  Define  storage  vectors. 
kmax  =  Tf/dt+  1; 
xm  =  zeros(4,kmax); 
xt  =  zeros(4,kmax); 
time  =  zeros(l  Jcmax); 
sig  =  zeros(3daiiax); 
am  =  zeros(l,kmax); 
um  =  zeros(2,kmax); 

Rtm  =  zeros(l,kmax); 
xm(;,l)  =  [0;0;0;2000]; 
xt(:,l)  =  [1000;0;10000;-3000]; 
Rtm(l)  =  sqrt((xt(l,l)-xm(l,l))^2  + 
acm  =  0.0; 
act  =  0.0; 


%  xm  =  [  y  yd  X  xd  ]' 
%  xt  =  [  y  yd  X  xd  ]' 


%  Initial  Conditions  for  Missile. 

%  Initial  Conditions  for  Target. 
',l)-xm(3,l))^2);  %  First  value  for  Range. 
%  Initial  acceleration  for  the  missile. 

%  Initial  acceleration  for  the  Target. 


%  Begin  the  Simulation  Loop, 
for  i  =  l:kmax-l 
%  Define  angles. 

beta  =  atan2(xt(2,i),  xm(4,i));  %  Velocity  angle  for  the  Target, 

gam  =  atan2(xm(2,i),  xm(4,i));  %  Velocity  angle  for  the  Missile. 

sig(l,i)  =  atan2(xt(l,i)-xm(l,i),  xt(3,i)-xm(3,i));  %  LOS  angle. 

sig(2,i)  =  ((xt(3,i)-xm(3,i))*(xt(2,i)-xm(2,i))  -  (xt(l,i)-xm(l,i))... 
*(xt(4,i)-xm(4,i)))/Rtm(i)^2; 

sig(3,i)  =  ((xt(3,i)-xm(3,i))*(act*cos(beta)-acm*cos(gam))-... 
(xt(l,i)-xm(l,i))*(act*sin(beta)+acm*sin(gam))+... 
2*(xt(4,i)-xm(4,i))*(xt(2,i)-xm(2,i)))/Rtm(i)^2; 
sig(3,i)  =  sig(3,i)+2*(xt(4,i)-xm(4,i)+xt(2,i)-xm(2,i))*sig(2,i)/Rtm(i); 


%  Acceleration  —  Max  value  perpendicular  to  gamma. 
am(i)  =  N*sign(sig(2,i)  +  sig(3,i)*abs(sig(3,i))/{2*N)); 
um(;,i)  =  [am(i)*cos(gam);  -am(i)*sin(gam)]; 
xm(:,i+l)  =  Phim*xm(;,i)  +  Delm*um(:,i); 
xt(:,i+l)  =  Phit*xt(:,i)  +  Delt*[0;0]; 
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tinie(i+l)  =  time(i)  +  dt; 

Rtin(i+1)  =  sqrt((xt(l,i+l)-xm(l,i+l))'^2  +  (xt(3,i+l)-xm(3,i+l))'^2); 
acm  =  am(i); 
end 

%  Evaluation  of  Miss  distance, 
iflag  =  0; 

rm  =  find(Rtm==inin(Rtm)); 
if  rm  ==  kmax 
iflag  =  1; 

It  =  rm  - 1; 

elseif  Rtm(rm+ 1  )<Rtm(rm- 1 ) 

It  =  nn; 
else 

It  =  rm  -  1 ; 
end 

%  Average  Velocities  in  Intercept  Area. 

Vm  =  .5*(sqrt(xm(2,It)^2+xm(4,It)''2)  +  sqrt(xm(2,It+l)'^2+xm(4,It+l)'^2)); 

Vt  =  .5*(sqrt(xt(2,It)^2+xt(4,It)''2)  +  sqrt(xt(2,It+l)^2+xt(4,It+l)^2)); 

QE  =  .5*dt*(Vm  +  Vt);  %  Quantization  Error, 

if  iflag  =  0 

r  =  interp(xm(:,It;It+l),xt(:,It:It+l));  %  Interpolated  minimum  miss  distance, 
rstr  =  ['Inter  Miss  =  '  num2str(r) '  ft']; 
else 

rstr  =  [’Crossover  never  reached']; 
end 


%  Index  of  the  minimum  range  value. 

%  Determine  whether  the  crossover  occurs 
%  before  or  after  the  min  range  value. 


%  Plot  the  output  of  the  simulation  with  results. 

plot(xm(3,:),xm(l,:),'-w',xt(3,:),xt(l,:),'-w’) 

text(.  15, .85, ['Intercept  Time  =  ’  num2str(time(rm)) '  sec’],'sc’); 

text(.15,.81,[’Min.  Miss  =  '  num2str(Rtm(rm)) '  ft’],’sc'); 

text(.15,.77,['Quan.  Err  =  '  num2str(QE)  ’  ft’],'sc’); 

text(.15,.73,rstr,'sc'); 
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text(.15,.69,['N  =  ',num2str(N)],’sc'); 
titleCControl:  2nd  Order  Switching  Laws  (sig,  sigd)') 
xlabeK'Direction  1  (ft)'),ylabel('Direction  2  (ft)') 
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4.  SIMPRNV.M 


%  SIMPRNV.M  simulation  of  the  missiWtarget  simulation. 

%  Proportional  Navigation  Guidance. 

%  written  by  Colin  R.  Cooper 

%  Define  states. 

As  =  [0  1;  -64  -16]; 

Bs  =  [0;  64]; 

Am  =  [0  1  0  0;0  0  0  0;0  0  0  1;0  0  0  0]; 

Bm  =  [0  0;10;0  0;0  1]; 

At  =  [0  1  0  0;0  0  0  0;0  0  0  1;0  0  0  0]; 

Bt  =  [0  0;l  0;0 0;0  1]; 

Tf  =  2.25; 
dt  =  .01; 
maxac  =  1000.0; 

%  Descretize  the  states. 

[Phis, Dels]  =  c2d(As,Bs,dt); 

[Phim.Delm]  =  c2d(Am,Bm,dt); 

[Phit,Delt]  =  c2d(At,Bt,dt); 

%  Create  strage  vectors, 
kmax  =Tf/dt  +  1; 
xm  =  zeros(4,kmax);  %  xm  =  [  y  yd  x  xd  ]' 

xt  =  zeros(4,kmax);  %  xt  =  [  y  yd  x  xd  ]' 

um=zeros(2,kmax); 
time  =  zeros(l,kmax); 

Rtm  =  zeros(l,kmax); 
b  =  zeros(2,kmax); 


%  Estimator  System. 
%  Missile  System. 

%  Target  System. 


%  Estimator  for  sigma-dot. 

%  Missile  State  Equations. 

%  Target  State  Equations. 

%  Simulation  Time. 

%  Sample  step  size. 
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xm(:,l)  =  [0;0;0;2000];  %  Missile  Initial  Conditions. 

xt(:,l)  =  [10(X);0;l(XXX);-3000];  %  Target  Initial  Conditions. 

Rtm(l)  =  sqrt((xt(l,l)-xm(l,l))''2  +  (xt(3,l)-xm(3,l))''2);  %  First  value  for  Range. 


%  Begin  simulation  loop, 
for  i  =  l:kmax-l 

Vm  =  sqrt(xm(2,i)''2  +  xm(4,i)''2); 

ub(i)  =  atan2(xt(l,i)-xm(l,i),  xt(3,i)-xm(3,i)); 

b(:,i+l)  =  Phis*b(:,i)  +  Dels*ub(i); 

gam  =  atan2(xm(2,i),  xm(4,i)); 

if  Vm*b(2,i)<=maxac  %  Apply  Tlirust  limitation. 

am  =  Vm*b(2,i); 
else 

am  =  maxac; 
end 

um(:,i)  =  [am*cos(gam):  -am*sin(gam)];  %  Acceleration  is  perpendicular  to  gamma. 
xm(;,i+l)  =  Phim*xm(:,i)  +  Delm*um(:,i); 
xt(:,i+l)  =  Phit*xt(:,i)  +  Delt*[0;0]; 
time(i+l)  =  time(i)  +  dt; 

Rtm(i4l)  =  sqrt((xt(l,i+l)-xm(l,i+l))'^2  +  (xt(3,i+l)-xm(3,i+l))'^2); 
end 

rm  =  find(Rtm==min(Rtm));  %  Index  of  minimum  range  value, 

if  Rtm(rm+ 1  )<Rtm(rm- 1 )  %  Determin  whether  crossover  occurs 

It  =  rm;  %  before  or  after  the  min  range  value, 

else 

It  =  rm  - 1 ; 
end 

r  =  interp(xm(:,It;It+l),xt(:,It:It+l));  %  Obtain  the  Interpolated  min  miss  distance 

Vm  =  sqrt(xm(2,It)^2  +  xm(4,It)^2); 

Vt  =  sqrt(xt(2,It)^2  +  xt(4,It)''2); 

QE  =  .5*dt*sqrt(Vm^2  +  Vt''2);  %  Quantization  Error. 

%  Plot  the  output  of  the  simulation. 

axis([0  10000  0  1400]);  %  Same  scale  as  the  Optimal  Control  Sim. 
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plot(xm(3,:),xm(  1  ,;),'-w’,xt(3,:),xt(  1 

text(.  15, .85, ['Intercept  Time  =  '  num2str(time(rm))  ’  sec'],'sc') 
text(.15,.81,['Miss  distance  =  '  num2str(Rtm(rm)) '  ft'],' sc'); 
text(.15,.77,['Quant  Error  =  '  num2str(QE) '  ft’],'sc'); 
text(.15,.73,['Inter  Error  =  '  num2str(r)  '  ft'],'sc'); 
text(.15,.69,['Sigma(It)  =  '  num2str(ub(It)*  180/pi) '  deg'],'sc'); 
text(.15,.65,['Max  Acc.  Limit  =  '  num2str(maxac) '  ft/sec''2'],'sc'); 
titleCControl:  Proportional  Navigation') 
xlabeK'Direction  1  (ft)'),ylabel('Direction  2  (ft)'); 
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5.  BB3RDORD.M 


%  BB3RDORD.M  3/7/91 

%  This  is  a  simulation  of  a  3rd  order  system,  forward  tune. 

%  This  version  allows  varied  N  values. 

%  written  by  Colin  R.  Cooper 


xl0=-.5;  %  Define  initial  conditions 

x20=0; 

x30=0; 

N=1 ;  %  Set  the  magnitude  of  the  control  effort. 

Tf=2.5;  %  Set  maximum  time  of  simulation. 

dt=.002;  %  Set  the  simulation  step  size. 


A=[0  1  0;0  0  1;0  0  0];  %  Define  the  State  Equations 

B=[0  0  1]’; 

C=[l  0  0]; 


[PhiT)el]=c2d(A,B,dt); 
kmax=Tf/dt  +1; 

%  Prepare  storage  vectors. 

u=zeros(l,kmax); 

x=zeros(3,kmax); 

y=zeros(l,kmax); 

w=zeros(l,kmax); 

f=zeros(  1  ,kmax); 

time=zeros(  1  ,kmax); 


%  Discretize  the  system. 

%  Max  integer  value  for  the  simulation. 
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x(: ,  1  )=[x  1 0;x20;x30] ; 


%  Define  initial  conditions  in  state  vectors 


%  Begin  loop  for  simulation, 
for  (i=l:kmax-l) 

w(i)=sign(x(2,i)+x(3,i)*abs(x(3,i))/(2*N));  %  Defining  the  switching  law. 

f(i)=x(2,i)+w(i)*(x(3,ir2)/(2*N); 

u(l,i)=-N*sign(x(l,i)  +  (x(3,i)'^3)/(3*N'^2)  +  w(i)*x(3,i)*x(2,i)/N  +... 
f(i)*abs(f(i)/N)'^.5); 

x(;,i+l)  =  Phi*x(;,i)  +Del*u(l,i);  %  Calculate  the  state  values. 

y(l,i+l)  =  C*x(:,i+1); 

time(i- 1)=  time(i)  +  dt;  %  Store  time  vector, 

end; 

%  Plot  the  switching  law  and  its  components, 
clg,  axis([0  2.5  -1.5  1.5]); 

plot(time,u, time,. 75*w, time, f)  %  w  is  scaled  to  distinquish  it  from  u. 

title(’Control  Laws  vs  Time') 

pause 

%  Plot  the  3-Dimensional  view  of  the  simulation  from  45  deg.  azimuth 
%  and  45  deg,  elevation  angle.  (Pos.  xl  vector  is  out  of  the  screen 
%  and  towards  the  left,  Pos.  x2  is  out  of  the  screen  and  towards  right, 
clg,  axis; 

plot3d(x(l,:),x(2,:),x(3,:),45,45); 
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6.  ADJOINT.M 


%  ADJOINT.M  is  a  simulation  of  a  controlled  system  using  formard  time  simulation 
%  and  reverse  time  or  adjoint  solution. 


%  written  by  Colin  R.  Cooper 

A  =  [0  1  0:00  1;000]; 

B  =  [0  0  1]'; 

C  =  [10  0]; 
xlO  =  2; 
x20  =  0; 

N  =  1;  %  Maximum  control  effort. 

Tf  =  4; 
dt  =  .005; 

tsf  =  x20/N  +  sqrt(xl0/N  +  .5*(x20/N)''2); 
tsa  =  sqrt(xl0/N  +  .5*(x20/N)''2); 
kmax  =  Tf/dt  +  1; 

%  Define  the  Storage  Vectors 
X  =  zeros(3,kmax); 
y  =  zeros(l,kmax); 
time  =  zeros(l,kmax); 
imp  =  zeros(l,kmax); 
imp(l)  =  -N/dt; 

imp(round(tsf/dt)+l)  =  2*N/dt; 
x(l:2,l)  =  [xl0;x20]; 
[Phi,del]=c2d(A,B,dt); 

for  i  =  1  :kmax  -  1 

x(:,i+l)  =  Phi*x(:,i)  +  del*imp(i); 
y(l,i+l)  =  C*x(:,i+l); 


11  Mar.  1991 

%  Define  Forward  Time  State  Equations 

%  xl  initial  condition. 

%  x2  initial  condition. 

%  Length  of  simulation. 

%  Time  increment  for  simulation. 

%  Forward  Time  Switching  Time 
%  Adjoint  System  Switching  Time 
%  Max  length  of  storage  vectors. 

%  States  of  the  system 
%  Output  state  vector 

%  Impulse  vector  for  control  times. 

%  First  pulse  at  t  =  0 
%  Second  pulse  at  t  =  tsf 
%  Set  initial  conditions  for  x. 

%  Discretize  the  state  equations 

%  Begin  the  simulation  loop 
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time(i+l)  =  time(i)+dt; 


end 

miss  =  min(sqrt(x(l,:).''2  +  %  Find  the  min.  miss  distance 

tff  =  time(find(sqrt(x(l,:).^2  +  x(2,:).^2)==  ...  %  Find  the  time  of  the  min.  miss 

min(sqrt(x(l,:).''2  +  x(2,:).'^2)))); 

clg,plot(x(l,:),x(2,:).’-w’),grid  %  Plot  the  Phase  Plane 

title('Forward  Time  Phase  plot')  %  Lable  the  graph  and  display  the  desired 

xlabel(’Xl'),ylabel('X2’)  %  information 

text(.6,.80,['min  miss  =  ’  num2str(miss)],'sc’) 
text(.6,.77,[’tsf  =  '  num2str(tsf)  '  sec.’],'sc’) 
text(.6,.74,['tff  =  '  num2str(tf0  ’  sec.’],’sc’) 
pause 

%  Now  for  the  Adjoint  System, 
xa  =  zeros(3,kmax); 
time  =  zeros(l,kmax); 
impa  =  zeros(l,kmax); 
impa(l)  =  N/dt; 

impa(round(tsa/dt)+l)  =  -2*N/dt; 

[Phi,del]=c2d(A',C,dt); 

for  i  =  1  :kmax  -  1  %  Begin  the  simulation  loop 

xa(:,i+l)  =  Phi*xa(:,i)  +  del*impa(i); 
ya(l,i+l)  =  B'*xa(;,i+1); 
time(i+l)  ==  time(i)  +  dt; 
end; 

missa  =  min(sqrt((xa(3,:)-xlO).^2  +  (xa(2,;)-x20).''2));  %  Find  the  min.  miss  distance 

%  from  the  initial  conditions 

tfa  =  time(find(sqrt((xa(3,:)-xl0).'^2  +  (xa(2,:)-x20).''2)==  ...  %  Find  the  time  of  the 

min(sqrt((xa(3,:)-xlO).'^2  +  (xa(2,:)-x20).''2))));  %  min  miss  distance 


%  Define  storage  vectors 

%  First  impulse 
%  Second  impulse 
%  Discretize  the  Adjoint  System 


plot(xa(3,:),-xa(2,:),'-w'),grid 
tide('Adjoint  Solution  Phase  plot') 
xlabel('X3'),ylabel('X2') 
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%  Plot  the  Phase  Plane 
%  Label  the  graph  and  display  the  desired 
%  information. 


text(.6,.8,['min  miss  =  '  num2str(missa)],’sc') 
text(.6,.77,['tsa  =  '  num2str(tsa)  '  sec.'],'sc') 
text(.6,.74,['tfa  =  '  nuni2str(tfa)  ’  sec.'],’sc') 
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7.  ADJSIM.M 


%  ADJSIM.M  Adjoint  solution  to  missile  intercept  problem. 


%  written  by  Colin  R.  Cooper 

11  Mar.  1991 

n  =  4; 

%  Proportionality  Constant 

V  =  5000; 

%  Closing  Velocity 

B  =  [0  00  0  0  1]'; 

%  Define  State  Matrices 

C  =  [-100  0  10]; 

Tf=  1; 

%  Set  first  value  for  incremented  times 

dt  =  .01; 

%  This  outer  loop  runs  the  complete  forward  time  system  for  each  value  of  Tf,  storing 

%  the  output  into  vectors  at  the  end  of  the  loop. 

for j  =  1:4 

kmax  =  Tf/dt  +  1; 

%  Maximum  index  for  storage  vectors 

imp  =  zeros(l,kmax); 

%  Qeate  vector  for  impulsive  input 

imp(l)  =  1/dt; 

%  Define  the  pulse  at  t  =  0 

X  =  zeros(6,kmax); 

%  Create  storage  vectors  for  states 

y  =  zeros(l,kmax); 

time  =  zeros(l,kmax); 

%  Forward  Time  Simulation. 

count  =  0; 

%  Counter  to  indicate  the  computer  is  busy 

for  i  =  1  :kmax  -  1 

%  Begin  simulation  loop 

count  =  count  +  1 ; 

ki  =  l/(v*(Tf  -  time(i)  +  le-12)); 

%  1/Range 

A  =  [0  1  0  0  0  0 

%  Define  the  time  varying  A  matrix 

0  0  n*v  0  0  0 

000 100 
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-100*ki  0-100 -20  100*ki0 
00000  1 
00  0000]; 

Phi  =  eye(6)  +  A*dt  +  .5*A''2*dt'^2; 

x(:,i+l)  =  Phi*x(;,i)  +  B*dt*imp(i); 
y(:,i+l)  =  C*x(:,i+1); 
time(i+l)  =  time(i)  +  dt; 
if  count  ==  50 
disp('working'); 
count  =  0; 
end 
end 

TendO’)  =  Tf; 
yend(j)  =  y(length(y)); 
eval(['y',num2str(j),’  =  y;']); 
eval(['t',num2str(j),'  =  time;']);l 
Tf  =  Tf+  1 
end 

%  Adjoint  Simulation. 

Tf  =  4.25; 
kmax  =  Tf/dt  +  1 ; 
imp  =  zeros(l,kmax); 
imp(l)  =  I/dt; 
xa  =  zeros(6,kmax); 
ya  =  zeros(l,kmax); 
time  =  zeros(l,kmax); 
count  =  0; 
for  i  =  1  :kmax  -  1 
count  =  count  +  1 ; 
ki  =  l/(v*(time(i)  +  le-12)); 

A  =  [0  1  0  0  0  0 


%  Discretize  the  matrix  with  2nd  order 
%  expansion  -  drop  higher  terms 


%  Counter  indicates  computer  is  busy 


%  End  simulation  loop 
%  Record  final  time 

%  Record  corresponding  final  value  ymiss 
%  Save  output  vector 
%  Save  corresponding  time  vector 
%  Increment  to  next  Tf 
%  End  outer  forward  time  loop 

%  Set  a  simulation  time 
%  Maximum  index  for  vectors 
%  Create  vector  for  impulsive  input 
%  Define  the  impulse  at  t  =  0 
%  Create  storage  vectors 


%  Begin  simulation  loop 
%  1/Range 

%  Define  the  time  varying  A  matrix 
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0  0  n*v  0  00 
000100 

-100*ki  0-100 -20  100*ki0 
000001 
000000]; 

Phi  =  eye(6)  +  A'*dt  4-  .5*A''^2*dt'^2;  %  Discretize  the  matrix 

xa(:,i+l)  =  Phi*xa(:,i)  +  C*dt*imp(i); 
ya(:,i+l)  =  B'*xa(:,i+1); 
time(i+l)  =  time(i)  +  dt; 
if  count  ==  50 
disp('working'); 
count  =  0; 
end 

end  %  End  adjoint  simulation  loop 


plot(ti,yl,'-w',t2,y2,'-w',t3,y3,'-w',t4,y4,'-w',time,ya,’-w')  %  plot  output 

xlabeK'Time  (sec)’),ylabel('Ymiss  (ft)'), grid 
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8.  INTERP.M 


function  r  =  interp(xm,xt) 

%  INTERP.M  will  return  an  interpolated  value  for  the  min.miss  distance  given  the  states 
%  for  the  two  intercept  values. 

%  Written  by  Colin  R.  Cooper  26  Mar  1991 

%  Increase  the  sample  rate  by  a  factor  of  100  in  crossover  region. 

dax  =  (xm(3,2)-xm(3,l))/100; 

day  =  (xm(l,2)-xm(l,l))/100; 

dbx  =  (xt(3,2)-xt(3,l))/100; 

dby  =  (xt(l,2)-xt(l,l))/100; 

%  Define  the  storage  vectors  for  100  data  points, 
ax  =  zeros(l,100); 
ay  =  zeros(  1,100); 
bx  =  zeros(l,100); 
by  =  zeros(l,100); 

%  Set  initial  values  for  each  vector. 

ax(l)  =  xm(3,l); 

ay(l)  =  xm(l,l); 

bx(l)  =  xt(3,l); 

by(l)  =  xt(l,l); 

%  Assuming  a  straight  line  trajectory  with  constant  velocity  in  the  crossover  region,  create 
%  the  interpolation  data  sets, 
for  i  =  1:99 

ax(i+l)  =  ax(i)  +  dax; 
ay(i+l)  =  ay(i)  +  day; 
bx(i+l)  =  bx(i)  +  dbx; 
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9.  PLOT3D.M 


%  PLOT3D  is  a  3d  plotting  function  allowing  rotation  and  elevation  adjustments.  The  plot 
%  shows  a  3-D  curve  and  its  projection  onto  the  X-Y  Plane. 

%  X,Y,and  Z  data  must  be  passed,  and  if  no  azimuth  or  elevation  values  are  passed 
%  they  default  to  AZ  =  45°,  EL  =  30°.  The  Azimuth  is  the  angle  of  rotation  of  the  view 
%  angle  about  the  Z-Axis.  AZ  =  0°  is  looking  straight  down  the  X-Axis  at  the  Y-Z  Plane. 
%  The  elevation  is  the  angle  from  which  the  plot  is  viewed.  EL  =  90°  is  looking  down  the 
%  Z-Axis  at  the  X-Y  Plane.  The  elevation  angle  can  vary  from  -90°  to  90°. 

%  The  tick  marks  on  the  axis  will  default  to  10  marks  per  axis  and  the  values  will  be 
%  displayed  on  the  screen.  To  define  the  values  of  the  tick  marks  a  three  element  vector 
%  must  be  passed  containing  [dx  dy  dz]; 

%  Axis  values  will  be  calculated  and  fixed  unless  an  axis  vector  is  passed  to  the 
%  program;  [xmin  xmax  ymin  ymax]. 

%  The  Transformed  2-D  vectors  V  and  Vs  are  returned,  where  V  is  the  3-D  curve  and 
%  Vs  is  the  Projection  on  the  X-Y  Plane. 

% 

%  Example  :  [V,Vs]  =  plot3d(x,y,z,-45,30,dx,Ax) 

%  written  by  Colin  Cooper  4/26/91 

function  [V,Vs]=plot3d(x,y,z,az,el,dx,Ax) 

if  nargin  <  5,  el  =  30;  end 
if  nargin  <  4,  az  =  45;  end 


az=-az*pi/180;  el=el*pi/180; 

alpha  =  45*pi/180;  beta=30*pi/180;  %  Angles  for  mapping  onto  2D 

%  alpha  =  rot.  Beta  =  elv. 

Sl=[sin(az)  cos(az)  0 

-sin(el)*cos(az)  sin(el)*sin(az)  cos(el)]; 

S2=[sin(az)  cos(az)  0 

-sin(el)*cos(az)  sin(el)*sin(az)  0];  %  Trans  for  Projection. 
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[r,c]=size(x); 
ifr>  1 

X  =  x';  y  =  y';  z  =  z'; 
end 

V  =  Sl*[x;y;z]; 

Vs  =  S2*[x;y;z]; 

if  nargin  <  6 

dx(l)  =  (max(x)-inin(x))/10; 
dx(2)  =  (max(y)-min(y))/10; 
dx(3)  =  (max(z)-min(z))/10; 
end 

axl  =  [min(x)  max(x);  0  0;  0  0]; 
ax2  =  [0  0;  min(y)  max(y);  0  0]; 
ax3  =  [0  0;  0  0;  min(z)  max(z)]; 

axlt  =  [fliplr(0:-dx(l):min(x))  0:dx(l):max(x) 
zeros([0:-dx(l):min(x)  0;dx(l):max(x)]) 
zeros([0:-dx(l):min(x)  0:dx(l):max(x)])]; 
ax2t  =  [zeros([0:-dx(2);min(y)  0:dx(2):max(y)]) 
fliplr(0:-dx(2);min(y))  0:dx(2):max(y) 
zeros([0;-dx(2):niin(y)  0:dx(2):max(y)])]; 
ax3t  =  [zeros([0:-dx(3):min(z)  0;dx(3):max(z)]) 
zeros([0:-dx(3):min(z)  0:dx(3):max(z)l) 
fliplr(0:-dx(3):min(z))  0:dx(3):max(z)]; 

axl  =  Sl*axl; 
ax2  =  Sl*ax2; 
ax3  =  Sl*ax3; 
axlt  =  Sl*axlt; 
ax2t  =  Sl*ax2t; 
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ax3t  =  Sl*ax3t; 


minx  =  min([V(l,:)  Vs(l,:)  ax  1(1,:)  ax2(l,:)  ax3(l,:)]); 
if  minx  ==  0,  minx  =  -1;  end 

maxx  =  max([V(l,:)  Vs(l,:)  axl(l,:)  ax2(l,:)  ax3(l,:)]); 
if  maxx  ==  0,  maxx  =  1;  end 

miny  =  min([V(2,:)  Vs(2,:)  ax  1(2,:)  ax2(2,:)  ax3(2,:)]); 
if  miny  ==  0,  miny  =  -1;  end 

maxy  =  max([V(2,:)  Vs(2,:)  ax  1(2,:)  ax2(2,:)  ax3(2,:)]); 
if  maxy  ==  0,  maxy  =  1;  end 

clg 

axis(’square'); 
if  nargin  <  7, 

axis([minx-.3*abs(minx)  maxx+.3*abs(maxx)  miny-.3*abs(miny)  ... 
maxy+.3*abs(maxy)]); 
else 

axis(Ax); 

end 

b  =  axis; 
hold  on 

plot(axl(l,:),axl(2,:),'-w',axlt(l,:),axlt(2,:),'+w’,... 

ax2(l,:),ax2(2,:),'-w',ax2t(l,:),ax2t(2,:),'+w',... 

ax3(l,:),ax3(2,:),'-w',ax3t(l,:),ax3t(2,:),'+w') 

Iv  =  length(V); 

plot(V(l,:),V(2,:),'-w',Vs(l,:),Vs(2,:),'-w',... 

[V(l,l:15:lv);Vs(l,l:15:lv)],[V(2,l:15:lv);Vs(2,l:15:lv)],':w') 
plot([b(l)  b(2)  b(2)  b(l)  b(l)],[b(3)  b(3)  b(4)  b(4)  b(3)].’-w'); 
hold  off 

text(.22,.08,[’AZ  =  '  num2str(-az*  180/pi)  ’“’J.'sc'); 
text(.72..08,['EL  =  ’  num2str(el*  180/pi)  ’°’],’sc'); 
text(.22,.88,[’dx  =  '  num2str(dx(l))  ],’sc'); 
text(.50,.88,['dy  =  '  num2str(dx(2))  ],’sc’); 
text(.72,.88,['dz  =  ’  num2str(dx(3))  ],'sc’); 
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